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Abstract 

Several multigrid schemes arc considered for the numerical computation of viscous 
hypersonic flows. For each scheme, the basic solution algorithm employs upwind spatial 
discretization with explicit multistage time stepping. Two-level versions of the various 
multigrid algorithms are applied to the two-dimensional advection equation, and Fourier 
analysis is used to determine their damping properties. The capabilities of the multigrid 
methods are assessed by solving three different hypersonic flow problems. Some new 
multigrid schemes based on semicoarsening strategies are shown to be quite effective 
in relieving the stiffness caused by the high-aspect-ratio cells required to resolve high 
Reynolds number flows. These schemes exhibit good convergence lates foi Reynolds 
numbers up to 200 x 10 6 and Mach numbers up to 25. 
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1 Introduction 


Over the past few years the need for efficient numerical methods to solve the equations 
governing hypersonic viscous flows has become very obvious. Mostly, flow solvers used in 
current aerospace programs, such as X-30 or SAENGER, exhibit slow convergence towards 
the desired steady-state solutions, which leads to high computer costs, long turn-around 
times, and a slowdown in the efforts to design these vehicles. The reason for this is the 
appearance of flow phenomena with very different scales and with highly nonlinear behav- 
ior. We mention here the laminar and turbulent boundary layers at very high Reynolds 
numbers and their interactions with shocks and slip lines, and furthermore, shock/ shock 
interactions which generate complex flow fields. Many numerical techniques which were 
developed to assist convergence of subsonic or transonic flow calculations are found to be 
inappropriate for hypersonic flow applications. For example, the time step of many oxpli< it 
and implicit schemes, which allows the transient behavior of strongly nonlinear hypersonic 
flow phenomena to be captured, is highly limited. Consequently, thousands of time steps 
are needed to converge the thin boundary layers. 

A particular method which was successfully developed to accelerate convergence for 
a broad range of flow problems at both subsonic and transonic speeds is the multigrid 
approach. This method, which uses a sequence of successively coarser meshes in order to 
propagate disturbances throughout the flow field, combines nicely with explicit multistage 
time-stepping schemes [1]. Good convergence rates were obtained for inviscid flow and, 
later on, for viscous flows also [2-4], Initial attempts to apply this promising method to 
hypersonic flows failed for several reasons. Primarily, the shock capturing capabilities of 
the central- difference scheme used in [1-4] were found insufficient to resolve strong shocks. 
Subsequently, the shock detection mechanism built into the central-difference scheme was 
improved in [5,6], In order to have strong shocks and slip lines resolved with fewer com- 
putational points, the central-difference scheme was replaced with an upwind-type scheme 
in [7-9]- Since the high-frequency damping properties of upwind schemes are generally less 
controllable compared with central-difference schemes, a variant of the standard multigrid 
approach was also used in [8,9]. With this variant, additional coarse meshes are generated 
by semicoarsening in the different coordinate directions. This strategy was felt neccessaij 
to alleviate convergence problems associated with the high-aspect-ratio cells of the com- 
putational mesh. An additional problem encountered is that very high temperatures can 
occur in the stagnation region and near the surface at high Mach numbers; and hence, 
the time step of explicit schemes may be severely restricted by the viscous stability limit. 
It was found [6,8,9] that the viscous time-step limit can be removed by implicit residual 
averaging. 

It is worthwhile to notice additional published work on multigrid schemes for hy- 
personic flows. Decker and Turkel [10] analyzed the effect of boundary conditions and 
Runge-Kutta coefficients on multigrid convergence for hypersonic inviscid flows. Leclercq 
[11] analyzed two-level multigrid cycles with multistage schemes and upwind differencing 
in one dimension, and she presented two-dimensional computations of hypersonic flows on 
unstructured meshes. As a means to remove stiffness associated with high-aspect-ratio 
cells, Blazek et al. [12] introduced an upwind-biased form of the residual smoothing by 
which higher Courant numbers could be obtained. Thomas [13] used multigrid in combina- 
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tion with third-order upwind differencing and implicit approximate factorization schemes. 
Koren and Hemker [14] solved the steady Euler equations with multigrid and point relax- 
ation applied as the smoother. Using damped restriction and upwind prolongation, they 
reported impressive convergence rates for high-speed flows around blunt bodies. 

The present paper describes recent efforts to understand and to improve the use 
of multigrid schemes for the computation of hypersonic flows. First, various two- level 
niultigrid schemes with and without semicoarsening are introduced. Then we use Fourier 
analysis of the schemes, when applied to the two-dimensional convection equation, in order 
to study the behavior of their components. For each multigrid approach, the solver uses 
upwind discretization combined with an explicit multistage scheme. We next consider the 
numerical solution of the Navi or- Stokes equations for hypersonic flows. In Section 5, the 
basic elements of the flow solver for these equations are described. Some details concerning 
the application of the time-stepping scheme to fine and coarse grid problems are presented 
in the first part of Section 6. The extension of the two-level schemes to multilevel ones 
is then discussed. Elements of multigrid that are of particular importance for high-speed 
flow computations are given. In the results section, we consider three different hypersonic 
flow problems to assess the capabilities of the multigrid schemes. The effect of stiffness, 
arising from coordinate grids with high-aspect-ratio cells and from flow alignment, on 
the performance of the multigrid methods is examined. The benefits of semicoarsening 
are clearly demonstrated. Moreover, with the semicoarsening strategies being considered, 
good convergence rates are obtained for Reynolds numbers up to 200 x 10 6 and Mach 
numbers up to 25. 
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2 Multigrid Strategies 

To set the stage for the discussion relating to multigrid in subsequent sections of this 
paper, we first briefly describe the multigrid method and different execution strategies 
that will be considered. The multigrid approach is based on the Full Approximation 
Scheme of Brandt [15]. The grid transfer operators are those considered by Jameson 
[1]. Coarser meshes are obtained by eliminating alternate mesh points in each coordinate 
direction. Both the solution and the residuals are restricted from fine to coarse meshes. 
A forcing function is constructed so that the solution on a coarse mesh is driven by the 
residuals collected on the next finer mesh. The corrections obtained on the coarse mesh 
are interpolated back to the fine mesh. The multigrid schemes investigated within the 
present work are displayed in Figure 1. Figure 1 a) shows a two-level scheme with full 
coarsening. Restriction of the solution from the fine mesh, (m,n), to the coarse mesh, 
(m/2,n/2), is done by injection, whereas full weighting is used for the restriction of the 
residuals. Prolongation of the corrections is done by bilinear interpolation. Figure 1 b) 
shows a scheme with semicoarsening in the different coordinate directions. Again, injection 
and full weighting are used in the restriction process. The corrections obtained on the 
coarse meshes are averaged before adding them to the fine mesh. This is indicated by the 
numbers at the “up” arrows. Due to this averaging, half of the individual corrections on the 
coarse meshes is lost. It is, therefore, anticipated that the scheme in Figure 1 a) should be 
computationally more efficient, provided there is enough high-frequency damping obtained 
with the smoothing scheme of the fine mesh. In order to overcome this deficiency of the 
semicoarsening scheme, two more variants are considered. For the scheme of Figure 1 
c), the solutions on the coarse meshes are computed sequentially. Hence, the corrections 
obtained on the (m/2,n) mesh can be used to update the (m,n/2) mesh before time stepping 
(as indicated by the horizontal arrow). The sequential update of the second coarse mesh 
allows the full amount of corrections to be passed up to the fine mesh. Note that this 
multigrid variant is not compatible with the idea of parallel computations. An interesting 
compromise between schemes of Figures 1 b) and 1 c) was suggested by "Van Roscndale 
[16] (Figure 1 cl)). Here, only the corrections common to both of the coarse meshes, 
(m/2,n) and (m,n/2), are averaged, whereas the corrections to the modes living either on 
(m/2,n) or on (m,n/2) are passed to the fine mesh in full. This scheme does allow parallel 
computations for the coarse meshes. 


3 Fourier Analysis of the Scalar Advection Equation 

A crucial factor in constructing an effective multigrid method is the selection of a 
smoothing or driving scheme. Local mode (Fourier) analysis is generally applied to evaluate 
possible smoothers on the basis of stability and high-frequency damping properties. The 
screening of schemes is often performed with a smgle-gnd analysis. Since a stable single 
grid scheme may not be stable for the multigrid process, the behavior of a smoother with 
a particular multigrid strategy is needed. In addition, the multigrid process can have 
a substantial impact on the performance of the multigrid method. In fact, as we will 
demonstrate in this paper, semicoarsening can provide significant improvement, relative 
to full coarsening, in the damping of the multigrid, especially when there is a strong mesh 
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anisotropy due to high-aspect-ratio cells. 

A two-grid or two-level multigrid analysis has been applied by Jameson [1], Mulder 
[17], and Leclercq [11] using various schemes (i.e., multistage time stepping, different types 
of relaxation) for solving the Euler and Navicr-Stokes equations of fluid dynamics. In 
[1] Jameson introduces a multilevel uniform analysis and considers the linear advection 
equation in one space dimension. This approach represents a departure from the standard 
two-grid analysis given by Stuben and Trottenberg [18], which forms the basis for the 
analysis used in [17] and [11]. With the multilevel uniform analysis, fine-grid and coarsc- 
grid corrections are computed at all points of the fine grid. Then a nonlinear filter is 
applied to remove the coarse-grid corrections at fine-grid points not contained in the coarse 
grid. The filtering produces additional errors in the form of a carrier wave with a frequency 
depending on the fine-mesh spacing. This analysis does not allow for the coupling (aliasing) 
effects due to the restriction operator (fine to coarse grid transfer operator) in the multigrid 
method. However, it does offer the advantages of simplicity and application to more than 
two-level schemes. Thus, it allows the rapid comparison of multigrid algorithms. If a 
multigrid method is unstable or inefficient according to this analysis, then it is probably 
not a reasonable scheme. 

In this section, we consider the scalar two-dimensional advection equation. The mul- 
tilevel uniform analysis of [1] is extended to two space dimensions and applied to full 
coarsening and semicoarsening strategies. 

Consider an initial value problem governed by the scalar advection equation 

— w(x,ij,t) + Cw(x,y,t) = 0, (x , y) £ ft, (3.1) 

when' the linear operator is 


C = a 


d , d 

di + W 


a > 0, b > 0, 


the domain Q C 9ft 2 , and t £ 3ft + . Assume a periodic boundary condition for the scalar 
function ?e(r, y, f). Let fi = {(x,y) : 0 < x < 1, 0 < y < 1). Define a fine grid Gj and 
coarse grids G Cj i (l — 1,2) that cover the domain Q such that G c j C Gf. We generate 
grids G Ct i by eliminating every other mesh line of G / in one or both coordinate directions. 
First we describe the fine-grid discrete problem. Let the grid Gf contain mj x rij cells and 
have uniform spacings A xj and A yj. Let the discrete function reside at the Gj 

mesh point (iAx At each point (?, j)? wo consider a corresponding cell (Cj)ij with 
corners at (i-1/2, j -1/2), (z + 1/2, j-1/2), (* + 1/2, 7 + 1/2), and (i — 1/2, j + 1/2). Suppose 
we approximate the spatial derivatives of (3.1) with first-order upwind differencing. Then, 
we obtain 




t [(“’/>”> - («’/)?-.,/] - N, [(»;)", - . (3.2) 


where the Courant numhers are 


iVjc = A^cry, N,, = \ n (Jf 


4 


and 


= aAy/, X v = bAxf, aj = 


_ A t f 
AQ, f 


The superscript n denotes the time level n At. If we estimate the time step At / by 

NAQ f 


At f = 


a Ay f + bAx / ’ 


it follows that 


Nt = 


aN 

ci 4" bA / 


N„ — bNA ' 


a + b A f ’ 


(3.3) 


(3.4) 


where jV is the Courant-Frieclrichs-Lewy (CFL) number for an explicit time-stepping 
scheme, and A f is the cell aspect ratio A Xf/Ayf. Taking a Fourier transform of (3/2), we 
obtain 


A = - [N e f(0 ( ) + N„M,)] (»,)", 


(3.5) 


where 


f(8) = (1 — cos 8) -f- i sin 0, 
and i = y/—l. The transformed discrete function is 

m j — l 7i j — 1 

(^f)k u k 2 = J] 5/ ( w /)"i,>i ex pHM« + ii 0 v )] » 

h=o n=o 

where the phase angles 8 $ and 8 n are given by 

k, 


(3.6) 


8 c =27 r- 


m f 


8jj = 27 r 


n f 


with wave numbers 


1 


h = -(rm/ _ o m/ ’ k 2 = ” 1 )»" 

Equation (3.5) can be rewritten as 

M, = (w ■,)"+■ - (w/)” = -Z(« { , »,)(•»/)". 


;"/• 


(3.7) 


with Z being the Fourier symbol for the difference operator. Consider the explicit p-stage 
scheme 

(«7) (0) = ("/)" 

(u>/) (i) = (w f )° - a t a f (R f f- l \ l = 1,2, ■ • • ,p (3.8) 

(7C/)^ n + 1 ^ = 


5 



where Rf is a residual function defined as 


R/ = AQ jCfiVf. 

Using this scheme and (3.7), one can always represent the change in wj by 

Sibf = -F(9^,9 n )Z(0^6 v )(wf) n , 

and the symbol of the time-stepping operator F is given by 

F ( 0 £ , 0^ ) — Op c\ pOip—\Z + 2 Z 

— (— l) p (a>o>-ia>-2 - ■ - «i )Z p ~ l . 


(3.9) 


(3.10) 


(3.11) 


In (3.11) we have assumed that both the convective and dissipative parts of the symbol 
Z are evaluated on each stage of the time-stepping scheme. This scheme is a member of a 
general class of (p, q) schemes. The p refers to the number of stages, and q designates the 
number of evaluations of the dissipative contribution. Suppose we have a (p, 2) scheme. 
Let 

Zr — Re(Z ), Zi = i • Im(Z). 

Then, if p > 3, the F of (3.11) is replaced by 

F — a p — a p (a\Zn + a p _i Z i)Z^ + a p a p -i(aiZ r + o p _2 Zi)Z\ — ■ ■ • 

- ( — 1 y-\a v a p - x • • • a 3 )( ai Z R + a 2 Z r )Z^ 3 (3.12) 

+ (-iy- 1 (a p a p - 1 ---cn)ZZ r ’- 2 . 


In this paper, wc consider a (5,3) scheme where the dissipative evaluations are weighted. 
For this scheme, the symbol of the residual function corresponding to the ( k + l)th stage 
is written as 




<7 


-1 


Ziw (k) + z R 


1=0 


The weighting factors are as follows: 


E 

l=o 


Jkl 


1. 


(3.13) 


Too = 1 , 

7 io = 1 , Til = 0 , 

720 = T 3 , 721 = 0, 722 = 73 i 

730 =r 3 , 731 = 0 , 732 = 7 3 1 733 = 0, 

740 = r 3 r 5 , 741 = 0 , 742 = 73 r 5, 743 = 0, 744 = 7 5 i 


(3.14) 


where T 3 = (1 - 73), T 5 = (1 - 75)1 7 3 = 0.56, and -y 5 = 0.44. The symbol of the 
time-stepping operator is given by 


F = cv 5 [1 — o 4 Zi(l - « 3 Z/) — a i Z i Z 1 (a 3 Z 2 Zj — 7 3 Zr) — r 5 7 3 Z 3 Zfl] , (3.15) 
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where 


Zi=Zj + */ 5 Z R , 

— Zl + 1^Z R , 

Z?, = 02(1 — ociZi). 

One can extend the stability range of the explicit scheme of (3.8) with implicit residual 
smoothing. For two dimensions, the residual smoothing can be applied in the form 

(i - a*)(i - /?,v,a n ykf = n ( ’\ (3.ie) 


where the residual (7£/)^ is defined using (3.9) as 


(R/)"' = <7,(71,)™, 


(3.17) 


A and V are the usual forward and backward difference operators, and / = 1, • ■ • ,p, with p 
being the number of stages in the scheme. The variable coefficients (3 ^ and (3 tj are defined 


as follows: 


0$ = max 


/?, = max 


1 

4 

1 

4 


n i y 

N* 1 + 4' r n s ) 



N 1 
N* 1 + ^ r~£ 



(3.18) 


with N/N* the ratio of the Courant number for the smoothed scheme to that of the 
unsmoothed scheme, and r, ; c = A t) /A^ . A reasonable value for the parameter ip is 0.125. 
In practice, this implicit procedure allows the explicit stability limit to be increased by a 
factor of 2 to 3. For additional details concerning implicit residual smoothing, see [3,5]. 

If residual smoothing is applied on each stage of the time-stepping scheme, then the 
symbol Z(9^,9 n ) appearing in (3.10) and (3.11) is replaced with 


z{e^e n ) = S{ 9 ^ 9 n )z{ 0 ^e n \ 


S{9M = TJ x T-\ 


(3.19) 


where 

= i + 2/?{(i - cos^), r, = i + 2/?„(i - cos 9 , t ). 

Before considering the discrete problem on some coarse grid G c , i, we define the re- 
striction operators and their corresponding Fourier symbols. Assume that G c , i contains 
(m//2) x (nf/ 2) cells (full coarsening). Let T c j and Q c } denote the operators that transfer 
the fine-grid solution and residual to the coarse grid. Since any coarse-grid point belongs 
also to the fine grid, we have simple injection for the solution transfer operator. Thus, 

T c f ’ 1 wf = w f , T; 4 = 1. (3.20) 

To ensure the conservation property for the residual transfer, we define 


QfiRfkj = 


(3.21) 


7 


where f.i x and y y are the standard averaging operators for the x and y directions. The 
symbol of this operator is given by 


Qy 1 = (1 + cos^)(l + cos <?,,). (3.22) 

If we now apply the p-stage scheme , we obtain at the corresponding point of the G, i 
problem 


(« , c,i) (0) = (w f ) + = ( w f ) 

(u’c,i) (l) = ( w c,i)° ~ oiia c> i 

(w c , 1 ) (n+1) =K,i) (p) , 


+ _ 


(flcl)*'- 1 * + Pci 


1 = 1,2, — ,p 


(3.23) 


where the forcing function 


Pc,, = Q}(R,) + 


(3.24) 


and 

(i) f y = An ,Cf(w f ) + 
(n c ,,y°> = An c , l c c , I (,„ c , I )(°>. 

By successive substitution in (3.23), one can easily show that 



Sib c ,i = —Fc^Sc^o-c^Qy 1 (Rf) + 

(3.25) 

or 

&w c , l = —F Ct iS Ci ia Ct iQ C j:’ i crJ 1 Zfw'f 

(3.26) 

with 

(wf) + = ( wf) n -f 6wf. 

(3.27) 

Since 

1 

II 

<g 

(3.28) 

then 

(w f ) + = g f {w f ) n , 

(3.29) 

where 

the fine-grid amplification factor is 



9f — 1 ~ FfZf- 

(3.30) 


Moreover, in this case, where we are considering full coarsening, the fine-grid approximation 
is 

(ivj) n+l = ( iv f ) n + 6w f + 6w Ci i, 


and thus, 


( u >/) n+1 = 


1 — Fc,\S c ,\v c ,\Q C f l <? f 1 Zf 

9c9}{wf ) n , 


9f(w f ) n , 


(3.31) 
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with g c denoting the coarse- grid amplification factor. 

We now apply the multilevel uniform analysis to two semicoarsening multigrid strate- 
gies (see Figure 1). Let G c> 1 and G c , 2 be two coarse grids containing (m// 2) x n f and 
m f x (n//2) cells, respectively. In the case of semicoarsening with simple averaging, we 
express the Fourier transform of the update for the fine-grid solution as 

( Wf) n+1 — {wf) n + Siif + u)\8iv Ci i + w 2 8w C i2 , (3.32) 


where = u> 2 = 1/2. Then the amplification factor associated with the coarse grids is 
given by 

g c = 1 - ^ F c> iS c ,\<y c ,\Q C f l(T T Z f + ^, 2 *S'c, 2 ^c l 2 Q/ 2 */ *^/] . (3.33) 

With sequential semicoarsening, we can use weightings of 1.0 for the coarse-gnd corrections 
by improving the estimate of the initial solution on either grid G c>2 , as indicated in Figure 
1, or G c , i. The coarse-grid correction is given by (3.26), and the other one 8w c<2 is 

given by 


8w 


c, 2 


—F, 


c,2 


Sc, 2 <rc, 2 Q C f 2 ~ Z Ci2 F c ,iSc,i<yc t iQ c j l a f l Z f g f (w f) n . (3.34) 


By substituting (3.26) and (3.34) into (3.32), one can compute the amplification factor 

g(0£,8 v ). . 

Before we applied the two-dimensional multigrid analysis just discussed, we examined 

the stability properties for a large number of multistage time-stepping schemes, including 
the schemes published in [19,20]. We performed this preliminary study with the one- 
dimensional advection equation for a single grid. Figures 2 and 3 show two examples from 
this study. The dash line represents the locus of the symbol of the difference operator, 
and it must lie inside the absolute stability curve. The dissipative term for first-order 
upwind differencing is evaluated three times for both the three- and five-stage schemes. 
In the case of the (5,3) scheme, the dissipation is computed on the first, third, and fifth 
stages, and the weights are those given in (3.14). While the (3,3) scheme exhibits fairly 
good high-frequency damping, there is a substantial improvement in damping with the 
(5,3) scheme, at the expense of a little extra work. The coefficients for the (5,3) scheme 
are those determined by Tai [21]. 

In Figures 4 and 5, we present contours of the amplification factor g(8^,8ti) f° r the 
following cases: 1) single grid, 2) full coarsening, 3) semicoarsening with simple averaging, 
4) sequential semicoarsening. For each case, cell aspect ratios of one and ten are shown. 
The improvement in damping with the two-level schemes is evident. There is greater 
compression of g{9^8 n ) to the origin with the semicoarsening strategies. The aspect ratio 
(A) of ten contours indicate that the one- dimensional behavior of the driving scheme m 
the r/ direction is recovered for large A. They also show that modes associated with the / 
direction are clamped much better with the semicoarsening schemes. 

According to the current analysis, where the nodal point of interest is assumed to be 
common to all meshes being considered, the semi coarsening scheme with selective averaging 
has the same damping behavior as the one with simple averaging. In practice, however, 
we observe the expected improvement in damping when using selective rather than simple 
averaging. This will be demonstrated later. 
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4 Governing Flow Equations 

Let p,u,v,p,T,E , and H denote, respectively, the nondimensional values of density, 
velocity components in the x and y Cartesian directions, pressure, temperature, specific 
total internal energy, and specific total enthalpy. In addition, let e x and e y he unit vectors 
of the Cartesian coordinate system ( x , y), and let n be a unit vector normal to the surface 

5 enclosing a volume V. Then the two-dimensional, unsteady Navier-Stokes equations, 
neglecting body forces and heat sources, can be written in conservative form in a Cartesian 
coordinate system as 

~Q t + /,(* + • ndS = 0 , ( 4 . 1 ) 

where the solution vector W and the tensors E c , J- v are defined as 



' p 


pue x -f pve y 


0 


pu 


(pu 2 4 - p)e x + puve y 


® X “b 

w = 

pv 

T — 

5 J C — 

puve x + (pv 2 + p)e y 

i 

, w 
II 

T xy£ X + a y e y 


pE 


puHe x + pvHe y 


(ua x + vr xy - q x )e x + 

+(uT xy + vOy — q y )e y 


with 


(T T * *■' 


Gy — 


. du dv 
.'du dv 

A 1 a* + ai 


du 

~ 2 ^dx' 

dv 


T xy — Tyx 


/ du dv \ 
\5y + dx ) ’ 



# = £ + 

P 

The scaling factor Re 1 = ^ MRe~\ with M and Re representing the Mach and 
Reynolds numbers, respectively. In this paper, the working fluid is air, and it is assumed 
to be thermally and calorically perfect. That is, the equation of state is 


P = (7 - 1 )p(E - (u 2 + u 2 )/2), T = p/p 


(4.2) 
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The quantities p and A are the first and second coefficients of viscosity, respectively, and 
A is taken to be -f/i (Stokes hypothesis). Sutherland’s law is used to determine the 
molecular viscosity coefficient ji. The coefficient of thermal conductivity (fc) is evaluated 
using the constant Prandtl number assumption. The effect of turbulence is taken into 
account by using the eddy-viscosity hypothesis. In the present work, the turbulence model 
of Baldwin and Lomax [22] is used. 


5 Spatial Discretization 

The numerical approximation of (4.1) follows the method of lines, which decouples the 
discretization in space and time. The physical domain around the aerodynamic body is 
divided into quadrilateral cells by the generation of a body-fitted grid. The discrete values 
of the flow quantities are located at the vertices of the mesh cells. For the flux calculation, 
an auxiliary grid is used, which is defined by connecting the cell centers of the ongina ce s 
(see Figure 6). The integral equation (4.1) is approximated by the spatial discretization, 

yielding 1 


dt 


wu-^mj + cwu). 


(5.1) 


where V, denotes the area of the control volume surrounding the grid node (Q c),y 

and (Q t> ) represent approximations of the convective flux and viscous flux, respective y. 
The viscous fluxes are approximated by central differences using a local transformation 
from Cartesian coordinates to the curvilinear coordinates [4], In the following some de- 
tails of the upwind discretization of the convective terms are discussed. The mviscid flux 
through the interface, (i + 1/2, j), is evaluated as 


(Qc)j-j-i /2,j = 9 l(^c)«,j + c)*+l,j] ' S H-l/2,> + 2 ^ i + 1 / 2 b^‘+ 1 / 2 >J 


(5.2) 


Here S J+1/2 , is the surface vector of face (i + 1/2, j), and R is the right eigenvector 
matrix of the flux Jacobian in transformed space. Equation (5.2) separates the invisci 
numerical flux into the sum of an averaged term corresponding to central differencing and 
a dissipative term, which adapts the discretization stencil in accordance with local wave 
propagation. The flux function * is based on the second-order accurate upwind TVD 
scheme of Yee and Harten [23]. Here, we do not repeat the formulation for but simply 
indicate some important details in the present numerical evaluation of $. Second-order 
accuracy is obtained with the limiter 

__ Q, — l/2,j(Q»+l/2,/ + e ) + «i+l/2,j(Q'i-l/2,j 2 + g ) (5.3) 

9l ' 3 ~ Of j— i/2, j 2 + a »+l/2,j 2 + 2e 


where a is the first difference of the characteristic variable, R 1 AW, and e is a small 
constant to prevent division by zero. The flow quantities at the face (i + are 

evaluated by Roe’s averaging procedure [24], Note that the scheme is identical to Roe s 
first-order flux difference splitting for g = 0. Since Roe’s scheme may violate the entropy 
condition when the eigenvalues of the flux Jacobians vanish, the eigenvalues are modified 
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using an entropy function 'I'. This entropy correction has to be carefully designed for 
\ iscous flow calculations. The shear layers along solid walls are numerically smeared, if 
an entropy correction is applied to the eigenvalues associated with the convective waves. 
On the other hand, if cells with high aspect ratios are present, additional support for 
damping in the direction of the long side of a cell is needed in regions of low velocities, 
such as stagnation points. Therefore, the correction of the eigenvalues T is constructed as 
a function of the cell aspect ratio. For the convective waves, the function 


»(V) = 


and for the acoustic waves 


lVl, 2 ^ if|A^|>«5; 

+ (1 — /?)|A( 7 |, otherwise, 


(5.4) 


*(V) 


^ 2 ^~ S i otherwise. 


(5.5) 


In equations (5.4) and (5.5), A^ represents the / th eigenvalue of the transformed flux 
Jacobian in the ^-direction. The entropy parameter, 8, is given according to Muller [25] 

8 — <5A^(1 + (A,/A^) w ), (5.6) 

where A^, A ^ are the spectral radii of the flux Jacobians in the £ and rj directions, 0 < u> < 1, 

and 0.1 < 8 < 0.5. If U = ue x + ve y , and are directed areas associated with the £ 
and i? directions, and c is the speed of sound, then 

A< = IU • S* | + c|S € |, A, = |U • S„| + C |S,|. (5.7) 

The blending coefficient, (3, accounts for the cell aspect ratio. It is given as 

P = max(l - A^/A 7 ,k). (5.8) 

It is shown below that k should be zero for accurate computations of shear layers. Fur- 
thermore, we will demonstrate that a wide range of flow problems can be solved accurately 
with a single set of parameters, 6 = 0.25, u = 0.3, and k = 0. 

6 Time-Stepping Scheme for Hypersonic Flows 

The basic elements of the time-stepping scheme have been outlined in Section 3, 
and they are not repeated here. In the following sections emphasis is placed on recent 
improvements to enhance robustness for hypersonic flow calculations. 


6.1 Multistage Scheme for the Fine and Coarse Meshes 

Consistent with the results for the advection equation in Section 3, we have observed 
the need to pair spatial discretization and particular time-stepping schemes for the solution 
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of the Navier-Stokes equation. The most robust choice of spatial discretizations found to 
this point is to use the second-order upwind scheme of Section 5 on the fine meshes and to 
set the limiter to zero everywhere on the coarse meshes. An alternative choice taken in [6,8] 
is to use scalar second- difference dissipation terms on the coarse meshes. This turned out to 
be less robust because the second differences are less diffusive with respect to the acoustic 
modes; also, the central-difference scheme allows waves to travel upstream in supersonic 
flow. As indicated previously, a five-stage scheme with three evaluations of dissipation and 
the coefficients of Tai [21] is used for time advancement. The large viscous stability limit 
of this scheme is shown in Figure 3. Disturbances are most effectively expelled out of the 
computational domain by using local time stepping and implicit residual smoothing [5,6]. 
The smoothing of the residuals allows a Courant number ratio (C F L / C F L*) as high as 
2.5 when CFL* = 2.3, which is roughly the stability limit of the explicit scheme. The 
time step is determined by the spectral radii of the inviscid flux Jacobi ans in the different 
coordinate directions, and A, ; , as 

A t = CFL-^~. (6.1) 

A^ + A, ? 

In order to stabilize the schemes in regions where the viscous stability limit is more 
restrictive than the inviscid limit, the coefficients of the implicit residual smoothing oper- 
ator are locally increased, as outlined in [6,8]. At strong shocks, however, high Courant 
numbers are not appropriate. Consequently, an adaptive time step is employed. By using 
the nondimensional second difference of the pressure as a switch, the value of CFL is 
locally reduced to about two at the shock. 

6.2 Multigrid Schemes 

For the numerical solution of the Navier-Stokes equations, the two-level strategies 
presented in Figure 1 are extended to multilevel schemes, as displayed in Figure 7. The 
only differences between the two-level schemes and the multilevel schemes occur in the 
restriction process. Whenever two “down” arrows meet at a coarse mesh, averaging is 
used to obtain the restricted variable. The multilevel arrangement of the coarse meshes, 
shown in Figure 7 b), was first given by Mulder [17], who used semicoarsening in order 
to solve the flow alignment problem. Suitable coordinate meshes for thin boundary layers 
exhibit mostly cells with high aspect ratios in the surface-aligned direction. Figure 8 
displays further variants of semicoarsening for these situations which are computationally 
cheaper than the semicoarsening schemes shown in Figure 7. 

One may notice that the central restriction and prolongation operators discussed in 
Section 2 allow for upstream propagation of disturbances in supersonic flow. Further- 
more, the corrections given by the standard multigrid scheme near strong shocks lead to 
divergence of the calculation, particularly during the initial part of the transient phase. 
Therefore, the restriction operator is damped by using 

Rij = max(l — e, i / n \0)R,j, (6.2) 

where Ri t j is the standard restriction operator, and is a. switch to detect strong 

shocks, 

ei/ n) = max(ui,Vi + uVi-i,Uj,Vj + i,Vj-i), (6.3) 
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(6.4) 


"PiJ ~l~ Pi+l,j ^ Pi,j — 1 ^ pi f j 4~ 

Pi-1, j +2pi , j + Pi+1, 7 ’ J P»,i-1 +2p«,j +P.-.J+1 

The damping coefficient, k (n \ is given a value of about one in the start-up phase of the 
multigrid process and is decreased to a value of about 0.4 at later cycle numbers in order to 
allow for good asymptotic convergence rates. This is in line with the restriction damping of 
Koren and Hemker [14], who based their damping coefficients on a more physical analysis. 

A fixed V-type cycle with time stepping on the way down is used to execute the multi- 
grid strategies described above. The robustness of the overall scheme is much improved 
by smoothing the resultant coarse-mesh corrections before they are passed to the finest 
mesh. The smoothing reduces high-frequency oscillations introduced by the linear interpo- 
lation of the coarse-mesh corrections. The factored scheme equation (3.16) with constant 
coefficients around 0.1 is used for this smoothing. Also, the application of Full Multigrid 
(FMG) provides a well-conditioned starting solution for the finest mesh being considered. 

7 Numerical Results 

Three different hypersonic flow cases are used to assess the capabilities of the multigrid 
schemes. These are laminar Mach 10 ( M = 10) flow over a compression ramp, turbulent 
flow over a slender forebody at high Reynolds numbers, and laminar flow over an airfoil at 
high Mach number and high angle of attack. Table 1 gives a summary of the geometries 
and the flow parameters of the test cases. In this table, T , n j is the dimensional free- 
stream temperature, and T w is the specified wall temperature. Also, the finest grid used 
for each flow computation is characterized by the streamwise and normal leading-edge 
spacings (A s; e > A n/ e ), along with the normal spacing (A n te ) at the end or trailing edge of 
a geometry. 

The flow over the compression ramp is identical to Case 3.2 of the Workshop on 
Hypersonic Flows for Reentry Problems, Part II, held in Antibes, France, 1991. This 
allows comparisons with the performance of other computational methods published in 
[26]. Figure 9 displays the coordinate mesh generated for this test case. The low Reynolds 
number allows for a mesh with moderate aspect ratios between 5 and 50 near the wall. 
The 129 x 81 mesh is successively coarsened down to 9 x 6, which yields 9 grid levels with 
semicoarsening and 5 levels with full coarsening. It is expected that the semicoarsening 
strategy should eliminate most of the stiffness associated with aspect ratio . The converged 
flow solution is shown in Figure 10. The computed extent of separation in the corner is 
somewhat smaller for the coarse mesh than for the fine mesh. Note that the result of the 
fine mesh agrees very well with grid-converged computations published in [13]. 

In the next figures, we analyze the performance of the different multigrid schemes. 
For this purpose, computations are started from a solution which was converged to about 
plotting accuracy. Figure 11 compares the different schemes of Figure 7. The numbers 
indicate the final convergence rate (r) of the schemes and the rate of data processing 
(RDP) on a CRAY-YMP to advance one grid point by one multigrid cycle. It is seen that 
the sequential semicoarsening scheme (Figure 7 c)) gives by far the best convergence rate. 
For this scheme, the effect of the modifications shown in Figure 8 is investigated in Figure 
12. One finds that the meshes obtained by full coarsening and by semicoarsening in the 
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direction normal to the wall are both important to achieve good convergence rates. From 
Figures 11-12, we conclude that semicoarsening with a selected number of coarse meshes 
is most effective for this flow problem; however, full coarsening does a surprisingly good 
job because of its low work count. Figures 13-14 show improvements which may be gained 
by using more than a single time step on the coarse meshes. The full coarsening scheme, 
Figure 13, gives only marginal gains when using more than two time steps on the coarse 
meshes. The sequential semicoarsening scheme, Figure 14, gives an initially improved 
rate, whereas the final rate is not affected by more work done on the coarse meshes. It is 
thought that the capabilities of the multigrid approach are put to full use for this test case. 
Further improvements are forseen only if the remaining stiffness in the discrete equations, 
that is the differences in the characteristic speeds of acoustic and convective waves, can 
be overcome by some proper means. A comparison between single mesh and multigrid 
computations is given in Figure 15. We find that the multigrid scheme with sequential 
semicoarsening converges within one tenth of the computing time required for the single 
mesh scheme. 

The grid generated for the ramp flow is well suited to study the grid-alignment problem 
which occurs for inviscid flow over the ramp. Figure 16 shows convergence histories of 
various schemes obtained by using a slip- wall boundary condition and omitting the viscous 
terms in the governing equations. Generally, convergence is worse for the inviscid flow 
because the convective eigenvalues of the flux Jacobian in the normal direction are zero 
for most of the grid points. The second-order solution does not converge, regardless which 
multigrid strategy we use. With the flux limiter, equation (5.3), set to zero everywhere on 
the fine mesh, the solution converges, provided we use all the coarse meshes introduced by 
the semicoarsening approach. Note that this result is identical to the findings of Mulder 
[27]. We have to conclude that the present multigrid schemes cannot cope with the grid- 
alignment problem in the context of high-order spatial discretizations. 

The flow over a slender forebody is chosen to represent a generic configuration corre- 
sponding to a high-speed civil transport aircraft or an air-breathing space transportation 
system with low wave drag. The high Reynolds numbers yield thin boundary layers, which 
can only be resolved with highly clustered coordinate meshes and large aspect ratio cells. 
The mesh used for the present investigations is displayed in Figure 17. The cells near the 
wall have aspect ratios up to 25000. The flow computations were done with fixed transi- 
tion at 2 percent chord and with the assumption of an adiabatic wall. Figures 18-20 show 
the solution obtained on three successively refined meshes. Both the distributions of skin 
friction and wall temperature are accurately computed, even with just 25 points in the 
normal direction. The effect of numerical dissipation introduced by the entropy correction 
function equations (5.4)-(5.8), is shown in Figure 21. If the value of uj is increased to 
2/3, which is a typical value used in central-difference codes [2,4], the wall temperature is 
badly reproduced on the coarse mesh. Furthermore, if we introduce numerical dissipation 
for the convective waves in the direction of the short side of the cells by setting K > 0, 
the shear layers are numerically smeared; and hence, the wall temperature is adversely 
affected. We notice that a value of k = 0.1, which again is a representative dissipation 
level of current central-difference codes, yields reasonable solutions. However, the conver- 
gence behavior is not changed much with this relatively low level of numerical dissipation. 
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Hence, all the results presented for the three test cases were obtained with a single set of 
parameters, 6 — 0.25, u = 0.3, and k — 0. 

Next we investigate the convergence behavior of the multigrid schemes. The fine mesh 
with 257 x 97 points allows 11 grid levels to be used with semicoarsening. First we notice 
that it is not possible to run the full diamond-shape tree of coarse meshes. Obviously, the 
time-stepping scheme is not well suited to handle the extreme aspect ratios which occur 
on the 17 x 97 mesh, for example. Using the proper half of the diamond, which includes 
the meshes with relatively low aspect ratios, the numerical solution converges. Figure 22 
displays a comparison of the different multigrid strategies. The computations are started 
from a preconverged solution. The computer time to update a grid point by one multigrid 
cycle, RDP , and the final convergence rate, r, are also included. Again, the scheme with 
sequential semicoarsening converges best. The differences between the multigrid schemes 
for this case, having cells with very high aspect ratios, are larger than for the ramp flow. 
The final convergence rate of the scheme with sequential semicoarsening is 15 times better 
than the rate with full coarsening. A comparison of the performance for the complete 
FMG process is given in Figure 23. The sequential semicoarsening scheme takes 194 cycles 
and 570s CPU on a CRAY-YMP to reduce the averaged residuals to 10~ 2 on the fine 
mesh. The scheme with full coarsening takes 1024 cycles and 1430s, and the single mesh 
code takes 7762 time steps and 6190s to achieve the same convergence level. Note that 
residuals of 10“ 2 correspond to a solution which is converged within plotting accuracy. If 
we compared computer times to reach lower levels of residuals instead, the results would 
have been even better for the multigrid scheme with semicoarsening. 

Laminar flow over an airfoil at M = 25 and a — 30 degrees is chosen as a final 
test case in order to demonstrate that the multigrid method used here can handle very 
strong shock waves and highly expanded flow. Figure 24 shows the 257 x 81 mesh . The 
numerical flow solution, which is represented in Figures 25-28, features a large separated 
flow region with two distinct vortices. The resulting shear layers are not well aligned 
with the coordinate mesh; and hence, considerable numerical smearing is expected in 
those regions. The difficulties in resolving this highly separated flow are illustrated by a 
comparison of the results obtained from meshes with different grid densities. Obviously, 
the mesh with 129 x 41 points is still too coarse to resolve the separated flow region. The 
convergence history of the sequential semicoarsening scheme is shown in Figure 29. The 
residual drops 8 orders of magnitude within 400 multigrid cycles. Note that the largest 
values of dp/di are located at the shock wave. The solution converges to plotting accuracy 
within 100 multigrid cycles. 

8 Conclusions 

New multigrid schemes for hypersonic flow computations have been investigated. The 
basic solution algorithm employs upwind discretization and explicit multistage time step- 
ping. Various multigrid schemes with semicoarsening are introduced in order to overcome 
stiffness resulting from the mesh cells with high aspect ratio w'hich are necessary to resolve 
viscous flows. The basic components of the algorithm are examined with Fourier stabiliy 
analysis applied to the two-dimensional aclvection equation. Both the results of the Fourier 
analysis and the computations of high Reynolds number flows suggest that the semicoars- 
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ening approach is effective. For the first time, convergence rates for hypersonic viscous 
flows are shown which are similar or even better than those previously published for the 
transonic regime [3,4]. Further work is required to make the computational scheme less 
expensive. This is especially true for the coarse meshes used within the semicoarsening 
approach, which make up the major portion of the overall work count of the scheme. Fur- 
ther improvements of convergence rates seem possible if stiffness arising from the difference 
of characteristic speeds of acoustic and convective waves can be overcome. For this pur- 
pose, new techniques such as preconditioning of the flow equations [28,29] or characteristic 
residual smoothing [12] seem to hold promise. 
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c) Sequential semicoarsening d) Semicoarsening with selective averaging 
Figure 1: Two- level multigrid schemes investigated in the present work 
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(a) Stability curves with implicit residual smoothing (CFL = 2.6, /? = 0.75) 
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(b) Amplification factor with implicit residual smoothing (CFL = 2.6, j3 = 0.75) 


Figure 2 Stability plots for 3-stage Runge-Kutta scheme with 
first-order upwind approximation (coefficients : 0.105, 0.325, 1.0) 
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(a) Stability curves with implicit residual smoothing (CFL = 5.0, 8 = 1.0) 
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(b) Amplification factor with implicit residual smoothing (CFL = 5.0, (3 = 1.0) 

Figure 3 Stability plots for 5-stage Runge-Kutta scheme with first-order upwind approximation 
and 3 evaluations of dissipation (coefficients : 0.2742, 0.2067, 0.5020, 0.5142, 1.0) 


23 




(b) Two levels, full 


Figure 4 Contour plots of an 
scheme with first-order up 
of dissipation (coefficients 
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9 0 47 
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7 0.37 

6 0.31 
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(c) Two levels, semicoarsening with simple averaging, 
weights = 0.5, A = 1 (CFL = 5.0, CFL* = 2.4) 
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(d) Two levels, sequential semicoarsening, weights = 1.0, A = 1 (CFL = 5.0, CFL* - 2.4) 


Figure 4 Contour plots of amplification factor for 5-stage Runge-Kutta 
scheme with first-order upwind approximation and 3 evaluations 
of dissipation (coefficients : 0.2742, 0.2067, 0.5020, 0.5142, 1.0) 
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(a) One level, A = 10 (CFL = 5.0, CFL’ = 2.4) 
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(b) Two levels, full coarsening, A = 10 (CFL = 5.0, CFL = 2.4) 


Figure 5 Contour plots of amplification factor for 5-stage Runge-Kutta 
scheme with first-order upwind approximation and 3 evaluations 
of dissipation (coefficients : 0.2742, 0.2067, 0.5020, 0.5142, 1.0) 








(c) Two levels, semicoarsening with simple averaging, 
weights = 0.5, A = 10 (CFL = 5.0, CFL* = 2.4) 



(d) Two levels, sequential semicoarsening, weights = 1.0, A = 10 (CFL = 5.0, CFL* - 2.4) 


Figure 5 Contour plots of amplification factor for 5-stage Runge-Kutta 
scheme with first-order upwind approximation and 3 evaluations 
of dissipation (coefficients : 0.2742, 0.2067, 0.5020, 0.5142, 1.0) 
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c) Sequential semicoarsening d) Semicoarsening with selective averaging 


Figure 7: Multilevel schemes 


i 

i 


28 



29 





Figure 11: Influence of multigrid strategies on convergence for ramp-flow problem 



Figure 12: Influence of selected coarse meshes on convergence for ramp-flow problem 
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Figure 13: Influence of multiple time stepping on coarse meshes for multigrid with full 
coarsening 



Figure 14: Influence of multiple time stepping on coarse meshes for multigrid with sequen- 
tial semicoarsening 
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Figure 15: Convergence of ramp flow with single-mesh time stepping and multigrid with 
sequential semicoarsening 
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Figure 16: Convergence histories for flow- alignment problem 
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Figure 19: Distribution of skin friction along forebody 



Figure 20: Distribution of adiabatic wall temperature along forebody 
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a) Influence of the exponent u> 



b) Influence of the cut-off value, k , applied to the convective waves 

Figure 21: Influence of the entropy correction coefficients on the adiabatic wall temperature 
along forebody 
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Figure 22: Influence of multigrid strategy on convergence for forebody, mesh 256x96 



Figure 23: Convergence histories for single-mesh time stepping and multigrid with sequen- 
tial semicoarsening 
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Figure 24: Coordinate mesh for NACA 0012 
with 256x80 cells 



Figure 25: Mach contours for NACA 0012 
airfoil at M = 25, a = 30° 
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Figure 26: Streamlines in separated-flow region 
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Figure 27: Distribution of skin friction along airfoil 



Figure 28: Distribution of Stanton number along airfoil 
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